Budding and vesiculation induced by conical membrane inclusions 
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Conical inclusions in a lipid bilayer generate an overall spontaneous curvature of the membrane 
that depends on concentration and geometry of the inclusions. Examples are integral and attached 
membrane proteins, viruses, and lipid domains. We propose an analytical model to study budding 
and vesiculation of the lipid bilayer membrane, which is based on the membrane bending energy and 
the translational entropy of the inclusions. If the inclusions are placed on a membrane with similar 
curvature radius, their repulsive membrane-mediated interaction is screened. Therefore, for high 
inclusion density the inclusions aggregate, induce bud formation and finally vesiculation. Already 
with the bending energy alone our model allows the prediction of bud radii. However, in case 
the inclusions induce a single large vesicle to split into two smaller vesicles, bending energy alone 
predicts that the smaller vesicles have different sizes whereas the translational entropy favors the 
formation of equal-sized vesicles. Our results agree well with those of recent computer simulations. 

PACS numbers: 87.16.Dg, 87.17.-d, 82.70.Uv 



I. INTRODUCTION 

Cell membranes contain large amounts of proteins 
within or attached to the lipid bilayer The distribu- 
tion of the proteins is not necessarily homogeneous, which 
can have important functional consequences. For exam- 
ple, proteins with an intrinsic curvature couple to the 
bilayer conformation [5, 0, IH, d, 0, 13 5 ^^e one hand, 
such proteins are preferably found on similarly curved 
parts of the membrane Q, on the other hand, the pro- 
teins deform the membrane locally Asymmetric, 
curved proteins can regulate the polymerization of the 
three-dimensional cytoskeleton of the cell Il2l| an d con- 
trol intracellular transport via endocytosis [13l, [l3|. V irus 
endocytosis can occur via the same mechanism [l5|, [l6| . 
The conical inclusions in our model mimick asymmetric 
proteins within the bilayer j T,], p roteins or polymers at- 
tached to the bilayer 0, Ilq , curved lipid domains 
19| [201 , 1 211 1221] , and viruses that bind to the membrane 

The interaction between the inclusions in a lipid bilayer 
is mediated by membrane deformations and thermal un- 
dulations [23, 24], in addition to surface tension [2] and 
possible direct interactions that we do not consider in this 
paper. The deformation-induced, pairwise interaction of 
curved inclusions occurs in the absence of thermal mem- 
brane undulations and is usually repulsive |25L l26i] ; in a 
planar membrane it is long-range [2^, [2^, l27j|. However, 
the interactions can be strongly screened if the average 
curvature of the membrane and the protein curvature are 
similar |2^, [2^, [sO]- One obvious example for strongly 
screened interactions are inclusions that are placed on a 
vesicle with similar curvature radius [30]. Screening can 
also be achieved by many-body interaction in clusters of 
inclusions [2^, |3l|. At finite temperature, Casimir-like 
interactions due to membrane undulations generate at- 
traction [25, 26, 32, 33, 34, 35]. 

Curvature generation by inclusions and induced bud- 
ding in lipid bilayer membranes has been reported in 



many experimental studies of biological and biomimetic 
systems [3, [O, [IH, [3, . Computer simulations allow 
to study the membrane-mediated interaction between the 
inclusions in detail without the presence of other, direct 
interactions. Recently, bud formation by curved inclu- 
sions has been investigated with computer simulations 
|36|, l37|. It was found that the inclusions on the buds 
have a higher density than they had in the initially nearly 
fiat membrane [36]. This might appear to be a result of 
undulation-induced attraction that in consequence leads 
to clustering of the inclusions and to budding. 



Such systems and processes can be studied theoreti- 
cally on the basis of an elastic membrane that is char- 
acterized by its bending rigidity, and Gaussian sad- 
dle splay modulus, with curved inclusions that con- 
sist of sections of a sphere with a given opening angle. 
We demonstrate that bud formation can already be well 
understood on the basis of the membrane deformation 
alone. We show that the higher inclusion density on the 
bud is a result of a screened repulsive interaction. We 
further argue that the budding pathway plays an impor- 
tant role for bud size. This allows us to predict a range of 
possible bud radii for a given system, which nicely agrees 
with recent simulation results ^36i] . 



At finite temperature, the inclusions can exist in a fluid 
and in a crystalline phase, which depends on the strength 
of their repulsive interaction. We construct an approx- 
imate free-energy functional that takes into account for 
the bending energy as well as the translational entropy of 
the inclusions. We calculate a phase diagram for the fis- 
sion of a single vesicle of given size and for given number 
and geometry of the inclusions. The inclusion entropy 
plays a decisive role for the sizes of the smaller vesicles 
into which a larger vesicle may split. 
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II. MEMBRANE BENDING ENERGY 

A. Membrane shape near inclusions in a lipid 
bilayer 

The bending energy £^ of a lipid bilayer is given by the 
integral over the entire membrane area, 

S = J dS {2i^H^ + RK) , (1) 

where is the bending rigidity, R is the saddle-splay mod- 
ulus, H = {ci -\- C2) / 2 is the mean curvature, K = C1C2 is 
the Gaussian curvature, and ci and C2 are the principal 
curvatures at each point of the membrane. The integral 
over the Gaussian curvature is determined by the topol- 
ogy of the membrane and by the geodesic curvature at 
the boundary. In our case, the geodesic curvature is given 
by the geometry of the inclusions, so that in general this 
term of the integral over the membrane shape does not 
need to be calculated explicitly. For bud formation, we 
neglect the constant contribution of the Gaussian saddle 
splay modulus. 

In order to minimize the bending energy, the inclu- 
sions preferably order on a hexagonal lattice {Fig.^(a)); 
therefore it is a natural assumption that the symmetry 
axis is oriented normal to the local tangent plane of the 
vesicle on which the inclusions are placed. To calculate 
the deformation energy, we approximate the hexagons 
with overlapping circles that have the same projected 
area {Fig.\^(b)). 

If there are no overhangs, the membrane conformation 
can be described in Monge parametrization by a height 
field, hix^y)^ over a planar reference surface. For an al- 
most planar membrane, the bending energy of the mem- 
brane is 

£=\n j dA[^h{pf (p=(x,y)), (2) 

with J dA the integral over the reference plane. Min- 
imization of the bending energy gives the biharmonic 




FIG. 1: (a) Membrane deformations induced by curved inclu- 
sions in a planar membrane. The inclusions have a repulsive 
interaction potential that decreases with the distance between 
the inclusions, like V ^ d~^. To minimize the bending en- 
ergy, the inclusions order in an hexagonal structure, (b) The 
hexagons are approximated by overlapping circles that have 
the same projected area. 



Euler-Lagrange equation, 

A^h{p) = . (3) 

In cylindrical coordinates, the general solution of Eq. (|3]) 
is 

hip) = \p\2C2 - C3) + C4 + (Ci + \p^c^) Hp) (4) 

with the four integration constants Ci to C4 [38] . 

The boundary conditions that are imposed on the 
membrane are sketched in Fig. [2l The radius of the inner 
boundary, pi = Vi sin (a), and the slope of the membrane 
at the inner boundary, h'{pi) = a = — tan(a), are de- 
termined by the inclusion geometry. For n ~ AttR^cf 
inclusions on a vesicle with radius R and surface num- 
ber density a of the inclusions, the radius of the outer 
boundary is po ~ i?sin(/3) with jS = arccos((n — 2)/n); 
the slope of the membrane at the outer boundary is 
h'{po) = b = — tan(/3). For inclusions on a planar mem- 
brane, the latter expressions simplify to po ^ l/(7rcr)*^^/^^ 
and 6 = 0. The remaining two boundary conditions 
are given by fixing the membrane height at the inner 
(or equivalently at the outer) boundary and minimizing 
the energy with respect to the height of the inclusion 
above the vesicle (i. e. the height difference between both 
boundaries), which implies h{pi) = at the inclusion and 
dpAh{p)\p^ = at the outer boundary. 

Eq. dlj) together with the boundary conditions gives 
the shape of the deformation, 

7 / X _ (p^ - Pi){bpo - api) + 2poPi{apo - hpj) \n{p/pi) 

~ ^{Pi - P\) 

(5) 

and the corresponding bending-energy cost. 




FIG. 2: (Color online) Curved inclusion (red) and result- 
ing membrane deformation (blue). The inclusion geometry 
is characterized by the curvature radius, r^, the opening an- 
gle, a, and the projected inclusion radius, pi — ri sin(a). The 
size of the corresponding membrane patch is po, the slope 
of the membrane at the inclusion is a = — tan(a), and the 
slope of the membrane at the outer boundary is h (6 = for 
inclusions on a planar membrane). 
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The energy a function of po and 6, which depend on the 
inclusion density, while all other quantities are intrin- 
sic properties of membrane and inclusions. For a sin- 
gle inclusion in an infinite planar membrane, 6 = and 
Po CO, the bending energy vanishes and the membrane 
deformation is catenoid-like, h{p) = apiln{p/ pi). Note 
that in a pairwise approximation, the interaction energy 
for two inclusions in a planar membrane {b = 0) decays 
like d~'^ for large distances between the inclusions (large 
Po = d/2). 



B. Optimal, low, and high inclusion density 

For inclusions on a vesicle, the membrane shape and 
the minimal bending energy (assuming that the inclu- 
sions have maximal mutual distances) can be calculated 
using Eqs. ([5]) and (|6]). For bpo = api^ the mem- 
brane around the inclusion has almost catenoid shape 
[40]; the catenoid is a minimal surface without bending- 
energy cost. If the entire vesicle is covered with inclu- 
sions and catenoids such that the bending energy is zero 
(Fig. [3] (h))^ the inclusions have optimal density. 

For lower inclusion densities, in a first approximation 
the catenoid shape borders on a spherical shape with the 
curvature radius of the vesicle. The bending energy of 
a vesicle that is decorated with curved inclusions is re- 
duced by the fraction of the sphere's surface area that is 
covered by the inclusions and the catenoid-shaped mem- 
brane segments, see Fig. [3l Therefore for low inclusion 
density, the bending energy of the decorated spherical 
vesicle is in the range < 8 < Sttk,. In the full solu- 
tion, which is given by Eq. ([5|) and will be used in the 
remainder of the paper, there is no jump in the mean 
curvature from H = OtoH = l/R between the catenoid 
and a sphere as sketched in Fig. [3] but rather a smooth 
transition from zero to finite mean curvature. 




FIG. 3: (Color online) (a) Vesicle decorated with curved in- 
clusions. Around each inclusion, the membrane can be mod- 
eled by segments of the catenoid minimal surface (white). 
The total bending energy is ^ = S7Thi{l — Scat/Ssph), where 
Scat/Ssph is the area fraction of the vesicle that is covered 
with inclusions and catenoidal patches, (b) Vesicle decorated 
with inclusions at optimal density; the bending energy of the 
lipid bilayer membrane vanishes. 



For inclusion densities that are higher than the optimal 
density, due to the boundary conditions no solution can 
be constructed by matching of catenoids. In this case, 
the bending energy always has a finite value that can 
exceed the bending energy of a bare vesicle. 

The minimal bending energy of a vesicle with inclu- 
sions is shown in Fig. [4] as function of the number of 
inclusions, n, and the vesicle radius, R. We find degener- 
ate zero-energy ground states that have optimal inclusion 
density with an approximately linear dependence R{n), 
where ^] 

R^ \a\pin/4 ^ l/{7rcr\a\pi) = (cosa)/(7rcr sin^ a){l/ri) 
and (7) 
A/ (iKja^ pI) = (4cos^ a)/(7rcr sin^ a)(l/r^) 

The natural spontaneous curvature of the bilayer for 
given inclusion density and geometry is cq = ^ 
TiG\a\pi. 

The same value for the inclusion density can be high, 
optimal, or low, depending on the radius of the vesicle 
on which the inclusions are placed. The smaller the ra- 
dius of the vesicle, the larger the value of the optimal 
density. In a planar membrane, the slopes of two ad- 
jacent catenoid-like deformations cannot be matched for 
any finite distance between the inclusions. Therefore the 
inclusions are always in the high-density regime in this 
case. 




FIG. 4: Normalized bending energy, E/n^ of a vesicle with 
radius R with n inclusions {vi ~ 5.5 nm, a — 0.64). There is 
a region of low inclusion density at large R with < ^ < Sttk, 
which is delineated by a line of zero-energy ground states 
from a region of high inclusion density at small i?, where also 
bending energies 8 > Stvk, can be found. (The high energies 
that are cut off at small n and large R mark the breakdown 
of the small- curvature expansion of the bending energy.) 
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C. Budding and vesiculation 



n_) vesicles: 



Bud formation does not occur for a vesicle with low 
inclusion density and bending energy, < £ < Stta^, be- 
cause this would lead to an increase of the total bend- 
ing energy [42|. However, for high inclusion density 
(n ^ 4:R/{\a\pi)^ see Eq. ([7j)), the system can always 
reach a state of lower bending energy if small vesicles 
bud from the main vesicle. The set of smaller vesicles into 
which a large vesicle with high inclusion density splits up 
is not uniquely determined from bending energy alone, 
because the states of vanishing bending energy are de- 
generate. A natural assumption is that the vesicle will 
split into one large 'mother' vesicle and one or more small 
'daughter' vesicle(s) of equal size, such that the total 
bending energy vanishes and the membrane area is kept 
constant. In Fig. [5l we show the radii of the mother 
and daughter vesicles as function of the number of inclu- 
sions. For a given number of — 1 daughter vesicles, 
there is a maximal number of inclusions nmax = nl^'^wR 
{w = 4:/{\a\pi)) that still allows to obtain a zero energy 
state, for which mother and daughter vesicles have equal 
sizes. 

If the system can split up into riy smaller vesicles, it 
can also split up into a larger number of small vesicles 
[4^. For a vesicle with total number of inclusions n = 
n+ + {riy — l)n_ and radius R = {R\ + (n^ — 
bending energy minimization predicts for the radii and 
inclusion number on mother n+) and daughter (i^-. 




FIG. 5: (Color online) A single vesicle of radius R = 10 n 
with n inclusions splits into one large 'mother'-vesicle and 
several small 'daughter' vesicles, in the figure the cases of 1, 
2,3, and 4 daughter vesicles are shown. The same parameters 
as in Fig. |4] are used. For 100 inclusions, the initial vesicle 
has vanishing bending energy. For more than 100 inclusions, 
the sizes of mother and daughter vesicles are plotted. The 
upper branch always gives the radius of the mother vesicle, 
the lower branch is the size of the daughter vesicles. For a 
fixed number riv — 1 of daughter vesicles (as indicated), there 
is a maximum number of inclusions that allows the formation 
of a state with vanishing bending energy (for which mother 
and daughter vesicles have equal sizes; filled circles). 
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Because of the degeneracy of the states with vanishing 
bending energy, thermal fluctuations and the budding 
pathway play a decisive role to determine how a large 
vesicle with high inclusion density splits up into smaller 
vesicles. Note that the results of our analytical model so 
far do not depend on the value of the bending rigidity of 
the bilayer. 



D. Inclusion clusters 

A direct attractive interaction between inclusions can 
induce cluster formation [3]. In this case, the preferred 
curvature radius, Rq^ not only depends on the num- 
ber and geometry of the inclusions, but also on cluster 
size. For given inclusion curvature radius r^, opening an- 
gle Qfo, and fixed density ctq, inclusion clusters with a 
larger opening angle a and reduced density, a = (1 — 
cosao)/(l — cos a) ctq, have a stronger effect on the cur- 
vature of the membrane than homogeneously distributed 
inclusions. The preferred curvature radius for clusters 
with opening angle a is R{a) = ((1 — cos a) cos Q^)/((l — 
cos ao) sin^ a) / {TTViao) = f{a)/{7rriao). We call the nor- 
malized curvature radius f{a) the coagulation factor, be- 
cause it multiplies the curvature radius for the reference 




FIG. 6: The coagulation factor, f{a) = R(a)7Triao, describes 
the dependence of the optimal vesicle radius, i?o, on the de- 
gree of aggregation of the inclusions in the membrane. The 
total inclusion area is kept constant and the reference density 
(Jo corresponds to ao = 0.1 tt. 
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inclusion density and opening angle [45|,[4^, see Fig.[6l 
For a fixed number no of inclusions, the preferred 
curvature radius decreases when aggregates are formed, 
i. e. the efficiency with which the inclusions influence 
the membrane curvature increases. Cluster formation 
therefore also shifts the high-density regime to smaller 
numbers of inclusions for the same vesicle radius, no ^ 
f{a)AR/ri^ and may cause a large vesicle to break up 
into smaller vesicles. 



III. THERMAL FLUCTUATIONS 

For finite temperature, the translational entropy of the 
inclusions contributes to the free energy. We distinguish 
a crystalline hexagonal phase and a disordered fluid phase 
for which we construct free energy functionals. Phase di- 
agrams have been calculated more rigorously in Ref. [i^l 
in the limiting case of an almost planar membrane and 
for inclusions that are slightly stiffer than the membrane 
and weakly curved — but not in the context of budding. 

We neglect the interaction between inclusions by ther- 
mal membrane undulations. For a pair of inclusions and 
an arbitrary orientation of the inclusion axis, to lowest 
order of /dP the deformation energy is [25,, i2^, [l^, [3 

fdef. =8W^ (9) 

and the undulation-induced interaction energy [25|, [33 . 

^und. = -^^bT^ . (10) 

The ratio of the membrane deformation-induced re- 
pulsion to the undulation-induced attration in a planar 
membrane is Attko^ / {ZkBT). The undulation-induced 
attraction can be neglected if it is one order of magnitude 
smaller than the deformation interaction; for = 10 /c^T 
this is the case for a ^ 0.5, for k = 20 ksT already for 
a ^ 0.35. For low inclusion densities, the undulation- 
induced interaction energy decays with the square of the 
density while the free energy due to the inclusion entropy 
is of the order of ksT; in case of the phase diagrams in 
Fig. El the undulation free energy is about 10~^ ksT. 

A. Inclusion effective pair potential and effective 
hard-disc radius 

For inclusions on a hexagonal lattice, each inclusion 
corresponds to three pair interactions and a radius of 
the deformation patch, po, that is half the distance be- 
tween the inclusions. The effective pair potential, ob- 
tained from Eq. (|6]), is thus 

u{d) = '"ffr';^^^' if bd<2api ^^^^ 

u{d) ^ ' ifbd> 2api . 



For inclusions on a planar membrane {b = 0) and for 
large d^ i. e. for d ^ pi^ the repulsive interaction potential 
decays with a power law, u ~ d~'^. 

To determine the free energy of this system, we em- 
ploy the method developed for suspensions of repulsive 
colloids [iO], i- e. we mimic the interaction potential by 
hard discs with an effective radius, rhd- The radius is 
calculated from a comparison of the membrane deforma- 
tion energy with the thermal energy, ksT. We use a 
modified Barker-Henderson method that is appropriate 
for soft potentials [46, 50], 

ri,d = ^£" {l-exp[-/3S{p)]}dp, (12) 

where the upper integral boundary is determined by 
u{2pu) = ksT [51I, [52:]. The effective hard-disc radii 
therefore depend on the geometry of the inclusion, the 
bending rigidity of the membrane, and on the radius of 
the vesicle on which the inclusions are placed. 

In Fig. [7J the effective hard-disc radii, rhd, are plot- 
ted for inclusions with various opening angles as func- 
tion of the vesicle radius, R (an extremely large radius 

= 100 /im of the vesicle is used to describe inclusions 
in planar membranes). The hard disc radii increase with 
increasing vesicle radius; the increase of Vhd with R is 
the stronger, the larger the opening angle a of the inclu- 
sion is. For example, the inclusions with = 5.5 nm and 
a = 0.4 TT can approach each other by diffusion about 
an order of magnitude closer on a vesicle with radius 

= 10 nm than this is possible on a planar membrane. 
Consequently, the translational entropy of the inclusions 
lowers the free energy on the vesicle compared with the 
planar membrane. 

As discussed in Sec. IIIDl cluster formation of inclu- 
sions increases their effect on the membrane curvature. 




10 100 1000 10000 

R/r, 

FIG. 7: Effective hard-disc radii for inclusions on a vesicle as 
function of the vesicle's radius, R (see Eq. (|12|) ). All inclu- 
sions have the curvature radius, Vi — 5.5 nm, and the open- 
ing angles a = 0.47r, a = O.SStt, a = COStt (from top 
to bottom). The projected inclusion radii are in the range 
0.86 nm < p\ < 5.2 nm, the membrane bending rigidity is 
12 ksT. 
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a 

FIG. 8: Area fraction of effective hard discs for inclusions 
with n = 5.5 nm and opening angles 0.1 < a < 1.3. For a 
fixed number of inclusions, cluster formation leads to a larger 
area fraction of the effective hard discs and finally to crys- 
tallization. The area fraction of the effective hard discs with 
opening angle a is normalized by the area fraction of effective 
hard discs with a = 0.05 tt: an increase of a corresponds to a 
decrease of a (compare Fig. [6]). The inclusions are placed on 
vesicles with k = 12 ksT and various radii R = 0.1 /xm (dot- 
ted), R = 1 /im (short-dashed), R = 10 /xm (long- dashed), 
and R = 100 /im (solid). For comparison, the area frac- 
tion of effective hard discs is also shown for k — IksT and 
R — 100 /im (dashed-dotted). 



The effect of cluster formation on the area fraction of ef- 
fective hard discs is plotted in Fig. [HI For n = 12 ksT, 
which is a typical value for a lipid bilayer, clustering on 
vesicles with large radii, strongly increases the area 
fraction of the effective hard discs. To illustrate the 
strong effect of the bending rigidity on the effective hard- 
disc radius, which enters the calculation of the radius via 
the exponential function in Eq. ([T2|), we plot the effective 
hard disc radii for = IksT; the increase of the area 
fraction of the hard discs with a is much smaller than 
for = 12 ksT. Thus the translational entropy of the 
inclusions plays a much more important role for smaller 
[53,]. 



B. Inclusion entropy and free energy of hard discs 



gives the excess free energy. 



1 >^cs 

nksT 




- ln(l - y) . 



(14) 

The Carnahan-Starling excess free energy diverges at the 
crystallization transition of the effective hard discs. 

Because in the fluid as well as in the crystalline phase 
of the inclusions the squared thermal wavelength enters 
through the same constant and additive term, we consis- 
tently replace it in both cases — without loss of informa- 
tion — by the projected area of the inclusion, yrpf , such 
that J^id/n ^ kBT\og{(jTipl). 

Usually the translational entropy favors a homoge- 
neous distribution of particles. However, because the 
effective hard-disc radius depends on the membrane cur- 
vature, on a deformable membrane, a homogeneous dis- 
tribution of inclusion does not need to be the most 
favourable state. Instead, the inclusion density on the 
bud can be higher than on the mother vesicle because of 
the screened repulsive interactions. Fig. [9] shows the free 
energies of a fluid of effective hard discs with curvature 
radius = 5.5 nm for various opening angles in a lipid 
bilayer with n = 12 ksT. For nearly identical curvature 
radii of vesicle and inclusions, the effective hard-disc ra- 
dius almost coincides with the geometric hard-disc radius 
of the inclusion. 
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The free energy of a fluid of hard discs can be very 
well described by the Carnahan-Starling free energy [M, 
[HH. It is the sum of the ideal-gas free energy, Ti^jn ~ 
/i:BTlog(crA^), where A is the thermal wavelength (see 
e. g. Ref. [56]), and the excess free energy Tq^ [57]. The 
latter is calculated from the Carnahan-Starling equation 
of state ^H, 

with the hard-disc area fraction, y = cryrr^^. Integration 
of the thermodynamic relation p = —{OF / dV)T,N finally 



FIG. 9: Carnahan-Starling free energy of hard discs as func- 
tion of the inclusion area fraction. The system is depicted 
in the inset with periodic boundary conditions. We plot the 
free energy for discs with the geometrical projected radii of 
the inclusions, rhd = Pi = r^sina (solid line), as well as the 
free energies for effective hard discs for inclusions with curva- 
ture radius Vi = 5.5 nm and various opening angles on vesicles 
with K = 12 ksT: a = 0.16 and R = 20 nm (long dashed), 
a — 0.16 and R — 100 /xm (short dashed), a — 0.64 and 
R = 20 nm (long-dashed dotted), a = 0.64 and R = 100 //m 
(short-dashed dotted), a — 0.82 and R — 20 nm (dotted), and 
a = 0.82 and R = 100 /xm (double dashed). 
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C. Free energy per inclusion in fluid and 
crystalline phases 

We construct the free energy per inclusion in the crys- 
talline phase from the membrane bending energy and the 
fluctuation free energy of a harmonic crystal, and in the 
fluid phase from the sum of the membrane bending en- 
ergy and the translational entropy of the inclusions [58]. 

For the harmonic crystal, the spring constant /Cgp is 
obtained for a hexagonal lattice with the interaction po- 
tential in Eq. [TTl 

fep = fl^^TLs (3^^ + 4.p^i)pi {a^ -^b^)-{d^^ I2pl)ahd . 
[a ) 

(15) 

The free energy contribution of the positional fluctua- 
tions of the inclusions is therefore 



27r 



(16) 



or — after replacement of the thermal wavelength by the 
inclusion size — 



(17) 



Whenever we use the free energy for the crystalline phase 
in this paper, the Lindemann parameter remains below 
the critical value for melting of the inclusion crystal [60] . 

The transition between the fluid and the crystalline 
phase already occurs below the crystallization transition 
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FIG. 10: Free energies as function of the inclusion area 
fraction, crvrp^, for inclusions with n — 5.5 nm, a — 0.82 
in a membrane with k — 12 ksT for several vesicle radii: 
R — 1 jim. (solid), R — 100 nm (long-dashed), R — 50 nm 
(short-dashed), = 30 nm (dotted), = 20 nm (long dashed- 
dotted), R — 15 nm (short dashed-dotted). For low densities, 
the inclusions are in a fluid phase and the free energy is given 
by the membrane bending energy plus the Carnahan- Starling 
excess free energy for the effective hard discs. For high den- 
sities, the inclusions are in a crystalline phase and the free 
energy is given by the membrane bending energy and the free 
energy of a harmonic crystal. 



of the effective hard discs, cryrr^^ ^ 0.7. In Fig. [TOl 
the free energy per inclusion is plotted for several vesicle 
radii. Entropy reduces the optimal bud radius for a given 
inclusion density compared with Eq. ([7j). However, the 
bending energy alone still provides a good estimate for 
the optimal bud radius, because of the strong increase of 
the free energy per inclusion for low inclusion densities 
(see Fig.Iini). 



D. Vesiculation diagrams 



From the total free energy in Sec. IIII C[ we calculate 
vesiculation diagrams starting with a single vesicle of ra- 
dius R and a given number of inclusions for several val- 
ues of n. Because the topology changes when buds de- 
tach, the value of R plays an important role for vesicu- 
lation. Fig. [11] shows vesiculation diagrams for ^ = 
and for R = —hz/2 (the ratio R/k is still under debate, 
see Ref. [61]). We calculate whether the single vesicle is 
the most stable state or if fission in two or more smaller 
vesicles is favorable. With increasing initial vesicle ra- 
dius and inclusion density, fission becomes more likely to 
occur — first into two, at even larger R or a into three 
or more vesicles. 

For ^ = 0, fission in the bending-energy dominated 
regime produces two smaller vesicles that in general have 
different sizes, see Eq. ([8]). If entropy is important, the 
two vesicles may have equal size. In Fig. [11] (a), it is 
shown that a regime of equally-sized vesicles develops, 
bordering the three- vesicle regime and increasing in size 
with decreasing tz. Within the error bars of our calcula- 
tion, we find that the boundaries between one, two, and 
three vesicles are independent of the value of the bending 
rigidity for lOA^^T < < 200A:bT. 

For R = — /^/2, vesiculation takes place for smaller ini- 
tial vesicle radii and inclusion densities than for = 0, 
because each additional vesicle lowers the free energy 
by 4:7tR. In the entire two vesicle regime, both vesi- 
cles have different sizes for tz = 200 /c^T, = 100 /c^T, 
hi = bOksT, and = 30 ksT, while for = 10 ksT 
a small region of equally sized vesicles is found, see 
Fig. HI (b). 

Note that for bud formation, which has to occur be- 
fore vesiculation, the value of R is irrelevant and the phase 
diagrams for ^ = apply (assuming that the buds are 
connected by catenoidal necks with vanishing bending 
energy, and that the membrane area needed to form the 
neck is negligible). While the bud is being formed and 
has not yet detached, the integral over the Gaussian cur- 
vature and therefore the contribution of the saddle splay 
modulus to the deformation energy stays constant. How- 
ever, a negative saddle splay modulus facilitates the neck 
between two vesicles to break. In this case, the budded 
state can act as energy barrier for vesiculation that pre- 
vents fission, separating a high-energy single- vesicle state 
and a low-energy state of several smaller vesicles. 
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FIG. 11: (Color online) Vesiculation phase diagram for an 
inclusion density a of inclusions with ri — 5.5 nm and a — 
0.64 on (initially) a single vesicle of radius R. For small a 
or i?, the energetically favorable state is the single vesicle; 
fission, first into two, and finally into 3 or more vesicles is 
expected to occur when a and/or R is increased (more than 
3 vesicles are not resolved by the calculation), (a) R — ^: In 
the two-vesicle regime, a region where the two vesicles have 
equal sizes grows with decreasing n and bounds the three- 
vesicle regime. The lines depict the border between equally 
and differently sized vesicles for n — 200 ksT, k = 100 ksT, 
K, = 50/cbT, k = 30 ksT, k = lO/csT. (the vesicle sizes 
are not resolved in the 3+ region), (b) R = —k/2: Fission 
occurs already for smaller values of a and R. In the 2- vesicle 
region the two vesicles have different sizes for k = 200 /c^T, 
K = 100 ksT, At = 50 ksT, and k = 30 ksT. For k = 10 ksT, 
in a small region both vesicles have equal sizes. 



membrane, the initial membrane shape is decisive; for a 
fast membrane relaxation, the initial distribution of in- 
clusions mainly determines the budding process. 

The diffusion time is td = \^ where A is a char- 

acteristic length scale that the inclusion has diffused and 
D is the diffusion coefficient of the inclusion. A typical 
value \^ D = 1 /im^s~^ for the diffusion of lipids and the 
diffusion coefficient for proteins in cell membranes can 
be up to two orders of magnitude smaller The re- 
laxation time of the membrane is tr = rjX^ /{2tt^k) [g^I, 
where 77 = 1 mPa ~ 2.4 10~^^ kpT s nm~^ is the viscos- 
ity of the surrounding water and A is the wavelength of 
the membrane undulations. From the cubic versus the 
quadratic dependence on A, we find that for small A, 
tr <td^ while for large U > td- 

For hi = 10/c^T and D = 1 /im^s~^ (which is an up- 
per bound for the diffusion coefficient of proteins), we 
find that tr = td for A ^ 0.6 mm. This length is much 
larger than 10 /im, which is the size of cells [64] or giant 
unilamellar vesicles [65], thus the initial aggregation of 
inclusions is diffusion-limited. We therefore assume that 
inhomogeneities in the protein distribution on the mem- 
brane will immediately lead to a membrane deformation 
that minimizes the bending energy. The larger the initial 
inclusion density, the larger the spontaneous curvature of 
the membrane (compare Eq. ([7j)) and the smaller the size 
of the buds that are formed. 

Bud formation is initialized in some regions of the 
membrane that have a noticeably higher inclusion den- 
sity than others. The relative protein density fluctua- 
tions decrease with the size of a membrane patch which 
is considered. If we assume a random distribution of in- 
clusions, for a patch with N inclusions the relative fluctu- 
ations of the inclusion number are of the order of A^~^/^. 
Thus for small patches, the inhomogeneities are strongest 
and budding will preferably occur on the smallest possi- 
ble lengthscale. A small average membrane curvature 
with appropriate sign further attracts proteins to those 
regions where the bending energy per inclusion is already 
reduced. However, this clustering of inclusions during the 
budding process is hindered by a ring with opposite mem- 
brane curvature that forms the neck of a growing bud. 
This ring acts as energetic barrier that prevents further 
inclusions to enter a patch of the membrane where bud- 
ding has already started [1^ . Because of the neck forma- 
tion and the diffusion-limited budding process, the bud 
size is roughly determined by the initial inclusion density 
on the membrane. 



IV. BUDDING PATHWAY 



COMPARISON WITH SIMULATION 
RESULTS 



To shed more light on the role of the budding path- 
way, we compare the typical diffusion time of the inclu- 
sions with the relaxation time of the membrane confor- 
mation on the same length scale. If the diffusion of in- 
clusions is fast compared with the relaxation time of the 



Budding due to membrane inclusions has been studied 
recently with computer simulations [H, [s^ • In Ref. [s^l , 
entire vesicles with inclusions are simulated where the 
membrane is modelled as a dynamically triangulated sur- 
face. In Ref. [3S], coarse-grained model lipids are used to 
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FIG. 12: Energies per inclusion needed to place 10 inclusions 
of size pi — 2.5 nm on a vesicle with R—Vo nm, k — 2^ ksT, 
and R, = — 20/cbT, as function of a. Lines show bending 
energy (long- dashed), bending energy and inclusion entropy 
(short-dashed), saddle-splay energy (dotted), and total energy 
(solid). For comparison, we also plot the simulation data of 
Ref. [37,] (symbols indicate different calculation methods [37]), 
shifted by AF — —10.5 ksT (see main text). The deviation 
of the simulations and our theory for a ^ OA might be due 
to the surface tension used in the simulation, which is not 
included in our model. 



study budding for planar bilayer patches. 

In Fig. [121 the different contributions to the free energy 
per inclusion needed to graft 10 inclusions with given 
projected radius pi = 2.5 nm on a vesicle with radius 
R = 15 nm are plotted as function of the opening angle, 
a. The inclusions are in the fluid phase. For comparison, 
we also plot the simulation data of Ref. [ST^ , shifted by 
AF = — 10.5/c^T because our model does not account 
for thermal undulations of the membrane conformation. 
This energy difference is extracted from the simulation 
data for T = 300K and for T = 3K, see Fig. 6A in 
Ref. [37]. However, the very good match is somewhat 
fortuituous because we replace the thermal wavelength in 
the ideal gas free energy by the inclusion size, such that 
it is similar to the free energy obtained on a triangulated 
vesicle with a bond length that approximately equals to 
the inclusion diameter, compare [67|. 

We consider curvature radii of the inclusions that are 
both smaller and larger than the curvature radius of the 
vesicle. For a ^ and oo, such that pi = 2.5 nm, 

bending energy is needed to insert the flat inclusion into 
the curved vesicle. This bending energy cost decreases 
with increased opening angle a and is zero for a ^ 0.17 
where the curvature radius of the inclusion equals the 
curvature radius of the vesicle. If a is further increased, 
the bending energy per inclusion continues to decrease as 
more and more of the vesicle area consists of catenoidal 
patches around the inclusions. For a ~ 1.18, i. e. for 
even larger opening angles than plotted in the figure, the 
inclusions have optimal density and the bending energy 
gained by grafting all 10 inclusions to the vesicle is Sttk. 

In addition to the bending energy, there is an en- 



ergy cost for grafting that arises due to the saddle- 
splay modulus, which has been chosen to he R = —k, = 
—20 ksT for consistency with Ref. [s^]- The energy cost 
per inclusion only depends on the geodesic curvature of 
the membrane at the inclusion, i. e. on the opening angle 
a, which implies Sj^ = —27tR{1 — cos a). Therefore it is 
independent of vesicle radius and inclusion density and 
is only important to calculate the chemical potential for 
the inclusions on the surface; the budding transition at 
given inclusion density in the membrane is independent 
of R. 

In the simulations presented in Ref. [SG^], budding is 
induced by inclusions with = 5.5 nm that are initially 
placed in a regular array on a planar membrane with K: = 
12 ksT. Under the assumption that the initial inclusion 
density, a = 210~^nm~^, determines the bud size (see 
Sec. lIVp . we can roughly predict the bud radius from 
Eq. ([7j). Possible bud radii are estimated by comparing 
the free energies for different vesicle radii in Fig. [10] at 
the initial inclusion density. 

The parameters for which the free energies are plotted 
in Fig. [To] are chosen to match the bending rigidity and 
the inclusion geometry of the 'large inclusions' in Ref. [36*] 
with a = 0.26 TT [6S.]. For an initial inclusion density 
in the planar membrane, cr = 0.002 nm~^ ^ OAbairpf, 
which is estimated by visual inspection from the simu- 
lation snapshots, we find that the inclusions are in the 
crystalline phase on the vesicle with R = 1 pm (i. e. in a 
planar membrane). The free energy per inclusion is about 
10 ksT. Significantly smaller free energies per inclusion 
can be found for vesicle radii 22 nm ^ R ^ 100 nm; the 
optimal vesicle radii are 30 nm ^ ^ 60 nm with free en- 
ergies per inclusion of about —1 ksT. These radii agree 
well with the observed bud radius of R = 30 nm that 
formes in the simulations as final state via an initially 
slightly larger bud. The optimal bud radius from Eq. ([7j) 
is approximately 37 nm. 

Similarly, for the 'very large inclusions' in Ref. [3^ 
{a = 0.39 7r) we predict bud radii, 15 nm < < 20 nm, 
as observed in the simulations; based on bending energy 
only we find from Eq. ^ R ^ 11 nm. For the 'small 
inclusions' in Ref. [36] {a = 0.20 tt) that are already 
in the planar membrane almost in the fluid phase, our 
model predicts for the 36 inclusions studied in the simu- 
lations a maximal gain for the free energy per inclusion 
of ^ 1.5 ksT for R = 38 nm and a strong decrease of 
this energy gain for smaller vesicle radii. Already for 

= 34 nm and 29 inclusions, the energy per inclusion on 
the bud and in the plane are approximately equal. Larger 
energy gains are possible for larger bud radii, for which 
many more than the simulated 36 inclusions and a larger 
area of the bilayer patch are needed. From these consid- 
erations, it is not surprising that no budding is observed 
for the 'small inclusions' in the simulations of Ref. [36i] . 
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VI. CONCLUSIONS 

We have calculated the membrane-mediated interac- 
tion of conical inclusions in a lipid bilayer and the in- 
clusion entropy, which allow the prediction of budding 
transitions and vesiculation. Our model is based on the 
membrane bending energy; with this contribution alone, 
the spontaneous curvature of a bilayer with inclusions as 
well as budding can be predicted for many systems, rang- 
ing from protein inclusions to viral budding. Although 
the interaction between the inclusions by membrane de- 
formation is repulsive, the screening of the repulsive in- 
teraction due to the average membrane curvature allows 
higher inclusion densities on a bud than in the initial 
vesicle. Translational entropy of the inclusions favors the 
formation of equally-sized daughter vesicles and lifts the 
degeneracy that is found for states with vanishing bend- 
ing energy. 

From our calculations, the following picture of the 
effect of the bilayer deformation by curved inclusions 
emerges. For low inclusion density, the membrane around 
each inclusion assumes a catenoid shape of vanishing 
curvature energy. At optimal inclusion density, the 
catenoids are closely packed and the bending energy 
for the entire vesicle vanishes. For high inclusion den- 
sity, the boundary conditions for the membrane defor- 
mation around the inclusion do not allow the formation 
of catenoidal patches, and the inclusions always feel the 
membrane-mediated repulsive interaction with neighbor- 
ing inclusions. In this regime, bud formation can occur. 

If we assume a specific biological mechanism that leads 
to formation of clusters with well-defined and limited 
size, we find that such a mechanism can induce bud for- 
mation and vesiculation without the need to insert ad- 
ditional conical proteins into the cell membrane: clus- 
ter formation reduces the preferred curvature radius of 
the membrane. We quantify the effect of aggregation by 
the coagulation factor, which describes how the preferred 
curvature radius for a fixed amount of inclusions changes 
with the cluster size. 



In general, our analytical model is applicable for a wide 
range of length scales and inclusion geometries. Com- 
puter simulations are usually designed only for a specific 
length scale, e. g. a length scale comparable to the length 
scale of lipids in Ref. [36] or the lengthscale of entire vesi- 
cles in Ref. js^ • The good agreement with the simulation 
results of Refs. [36*, ^37] suggests that the approximations 
used in our calculations are justified. 

We argue that the undulation-induced attraction can 
be neglected compared with the deformation-induced re- 
pulsion and the translational entropy of the inclusions 
for f^a'^ > 15kBT/{27r), i. e. a> 0.35 for f<i = 20 ksT. 
For example, the BAR protein induces a local membrane 
curvature with an opening angle a ^ OA Jjl. Clathrin 
can induce a variety of opening angles [H, • 

The number of the inclusions per bud is determined 
by the budding process. Around a growing bud, a neck 
forms that presents an energetic barrier for the diffusion 
of inclusions. Because the deformation of the lipid mem- 
brane typically occurs much faster than the diffusion of 
the inclusions within the membrane, the number of inclu- 
sions per bud is well determined by the initial inclusion 
density in the membrane. From this, we can estimate a 
range of possible bud radii, which agrees well with the 
simulations in Ref. [36]. 

It would be interesting to test the validity of our model 
in the limits of (a) very small inclusions, such as lipids 
with large headgroups, when the description of the lipid 
membrane by curvature elastic constants may not be ap- 
propriate and (b) very fioppy membranes, when neglect- 
ing the thermal membrane undulations may not be jus- 
tified. 
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